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ABSTRACT 

We have studied the microscopic aspects of the spreading of liquid drops on a solid 
surface by molecular dynamics simulations of coexisting three-phase Lennard-Jones systems 
of liquid, vapor and solid. We consider both spherically symmetric atoms and chain-like 
molecules, and a range of interaction strengths. As the attraction between liquid and solid 
increases we observe a smooth transition in spreading regimes, from partial to complete to 
terraced wetting. In the terraced case, where distinct monomolecular layers spread with 
different velocities, the layers are ordered but not solid, with qualitative behavior resembling 
recent experimental findings, but with interesting differences in the spreading rate. 


INTRODUCTION 

The spreading of liquids on solid, surfaces has lately received considerable attention, due 
to its crucial role in materials processing and its ubiquity in zero-gravity environments, as 
well as the fascinating scientific issues which arise. Wetting has long been studied both theo- 
retically and experimentally [1-3], and many features are well understood within the context 
of the continuum equations augmented with van der Waals forces. However a number of 
significant microscopic issues are unresolved, such as the precise boundary condition appro- 
priate to a moving contact line and the dynamics of contact angle hysteresis, and for such 
problems molecular scale calculations can provide invaluable information. Previous work 
has examined the contact line problem [4], showing that the no-slip boundary condition 
breaks down in regions of unusually high stress, and several calculations in progress address 
the effects of heterogeneities and hysteresis, but here we will consider the somewhat related 
ill-understood phenomenon of “terraced wetting.” It has been observed, using interferomet- 
ric techniques [5], that certain polymeric liquids spread on atomically smooth substrates in 
the form of distinct monomolecular layers of nanometer thickness, which extend laterally 
over millimeter distances to completely wet the solid. Such films are too thin for a lubri- 
cation analysis to be applicable, and only heuristic alternative models based on Langevin 
dynamics or two-dimensional liquid layers are available. This situation motivated us [6] to 
conduct molecular dynamics (MD) simulations of spreading, based on the atomic scale mo- 
tion of the constituent atoms. While we are restricted to small drops, whose initial radius 
is only about 10 atomic sizes, our results complement the experiments by providing detailed 
three-dimensional information about the spreading dynamics. 
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SIMULATION METHOD 


In an MD computer simulation [7], one places atoms in a region of space with random 
initial velocities corresponding to the desired temperature, and integrates Newton’s equations 
of motion. The force on any one atom is due to interactions with the others, and we consider 
Lennard-Jones potentials of the form 


V(r) = 


0 



where the coefficients are related to the conventional Lennard-Jones energy and distance 
parameters e and a via a = <r -6 and /? = 4 ecr 12 , respectively. The fluid subsystem is used as 
reference, with all distances expressed in terms of <r = < 7 #, mass in terms of nif, and times in 
terms of r = (rrifcrg/eg) 1 ^ 2 . The parameter a = a&, the strength of the attraction between 
fluid and solid molecules, is varied from 1.0 to 1.4 to produce different wetting regimes, 
and Cff is kept at 1. For the solid-solid interactions we choose e„ = 50 and a m = 2~ 1 / 6 d, 
where d = \f2 is the lattice constant, so that the equilibrium position of the solid particles is 
also the minimum of the potential. The mass of the solid particles is 5mf. The potentials of 
fluid-fluid interactions and fluid-solid interactions are cut off at 2.5(7# and 2.5<7fo respectively, 
and the potential of solid-solid interactions is cut off at 

The “computational protocol” in our simulations is to consider 4000 fluid particles con- 
fined to a cube of side L = 60er, with the boundaries at z = 0 and z = L are replaced by 
a solid wall made of 5 layers of fee solid, totalling 9000 molecules, with the (100) surface 
exposed to the fluid. The solid is thick enough to prevent a direct interaction between the 
fluid particles on either side, and the remaining directions have periodic boundary conditions 
imposed. The equations of motion are integrated using the Beeman algorithm [7], with a 
time step 0.005r. The entire system is equilibrated at T = 0.7 with a = 0 and /?& = 1. 
The drop has then a roughly spherical shape of radius 10<7, surrounded by a cloud of vapor. 
The center of the drop is gradually moved to 16<r above the solid surface by time t — 35, 
whereupon the fluid-solid interactions attract the drop, leading to contact and subsequent 
spreading. The middle layer of the solid is kept at constant temperature T = 0.7 by rescaling 
of the velocities of the particles in that layer, in analogy to a laboratory experiment where 
the solid substrate is kept at a constant temperature. The computation ends when the drop 
attains a steady state, or reaches the edges of the simulation box. 


One difference betwen the present simulations and laboratory experiments is that the 
solid-liquid interaction cuts off at a fixed distance, whereas a realistic long-range van der 
Waals potential, arising from a superposition of many layers of r -6 interaction, behaves 
approximately as r -3 [8]. In order to check that the partial wetting is not due to this 
cutoff, we ran a simulation with an extended fluid-solid interaction. This extra potential 
had the form k/(z — z 0 ) 3 for z (the coordinate normal to the solid surface) greater than the 
cutoff, with the constants k and ro determined by the continuity of the potential there. We 
compared the final states with and without this extended potential for a = 1.0 and did not 
see any qualitative changes. 

We will also consider simple polyatomic molecules, made by introducing an attractive 
confining potential between adjoining atoms along a chain, e.g., Vi 2 (r) ~ r +12 + r +6 . If non- 
adjoining atoms retain a Lennard-Jones attraction, the resulting molecule is compact rather 
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than elongated, while appropriate repulsive or orientational potentials may be introduced to 
maintain a long chain. At the expense of further computational effort, other MD practitioners 
often work with realistic specific materials. 

MONATOMIC FLUIDS 

First we consider the spreading of a simple monatomic fluid, such as Argon, where the 
Lennard- Jones potential above is quatitatively valid. We studied a range of solid-liquid 
interaction strengths from a = 1.0 to 1.4, in steps of 0.1. For values a < 1.1 we observed 
partial wetting, in which the drop reached a stable shape, while for a > 1.2 the drop continues 
to spread up to the boundaries of the simulation box. In Fig. l(a,b,c) we illustrate the three 
types of behavior seen by showing the final states for a = 1.1, 1.2, and 1.4, respectively. 
The points in the figures are the centers of the fluid molecules, where the three-dimensional 
system has been projected onto the two-dimensional surface of the box. Note that the vapor 
density decreases as a increases, due to the condensation of vapor molecules onto the solid 
wall, and in Fig. 1(c) a condensed layer of regularly spaced molecules, commensurate with 
the solid structure, lies on the surface. There appear to be fewer molecules as a increases, 
but this illusion is caused by the enhanced ordering in the fluid. In Fig. 1(a) a true steady 
state is reached, with a well-defined contact angle. In contrast, in Fig. 1(b) the spreading is 
slow but does not terminate. Even at t — 885 the height of the drop is still decreasing while 
the particle “reservoir” in the bulk of the drop is nearly depleted - this case corresponds 
to complete wetting. In Fig. 1(c), the spreading is complete, and the drop evolves into two 
distinct molecular layers; we consider this case in detail. 

The time evolution of a terraced spreading drop is shown in Fig. 2, in the form of 
snapshots of the profile at various stages, with the final state given in Fig. 1(c). At t = 140, 
several liquid layers have formed, rather more prominently than in the a = 1.1 case, say, 
but layering in liquids near solid surfaces is not in itself unusual [8]. At t = 240 however, it 
is evident that the first two layers spread faster than the bulk of the drop, and that the first 
layer spreads faster than the second. Furthermore, a clear step-like structure has formed. 
At £ = 340 the spreading of the first layer continues, while the motion of the second layer 
has nearly come to a stop. At later times t = 440 and 540 the first layer continues to move 
outward while depleting the central bulk of the drop. 

The structures of the first and second layers in terraced spreading are quite strongly 
ordered. The particles form a defective hexagonal lattice, with stronger ordering in the 
first layer than the second, and within each layer the inner part is more ordered than the 
outer part. As the boundary of the layer moves outward, vacancies are created both in the 
interior and at the edge. When such plots are examined as a function of time, one sees 
that the vacancies in the first layer provide the likely sites for the particles to move in from 
above. Although ordered, the layers are by no means solid: since an MD simulations provides 
detailed information on the molecular motion, we have studied the intra-layer and inter-layer 
movement of particles in detail. Histograms of the horizontal and vertical displacements have 
a generally diffusive shape, with the trend that particles near the center of the first layer are 
most constrained, while the outer parts of the drop move most readily. In contrast to some 
theoretical models of terraced spreading, the conclusion here is that, except for perhaps the 
inner ring of the first layer, the liquid is not at all rigid or impenetrable. 
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For comparison with experiment, it is useful to quantify the rate of spreading of the drop. 
Fig. 3 shows the evolution of the average radii of the first and second layers as a function 
of time. Evidently, a reasonable fit is R 2 (t) = Clogt + D, with C « 430 and 272 for the 
first and second layers, respectively. Similar behavior is found for the other values of a 
where complete wetting occurs, with the same functional form and different values for the 
constants. This relation disagrees with the R 2 ~ t behavior found in laboratory experiments 
with non-volatile liquids [5] (“dry spreading”). One might suspect that the discrepancy is 
related to the fact that the condensation of vapor ahead of the drop reduces the degree of 
attraction of liquid to the surface. On the other hand, in this case the condensed vapor 
forms a distinct layer of fairly regular lattice, clearly separated from the spreading liquid, 
whereas only when a is smaller does the condensate mix with the drop. Another possible 
source of disagreement is the comparatively small number of molecules in the simulation. A 
laboratory drop has an enormous reservoir of molecules in its center available for continued 
spreading, whereas we have nearly exhausted the supply. Note however that the eventual 
leveling off of R 2 for the second layer resembles a simple finite-size effect, and suggests that 
a larger drop would continue to spread at the same rate. Indeed, a single run on a larger 
system, with 9000 fluid molecules and a maximum radius of 30, gives results consistent with 
those described above. 

CHAIN MOLECULES 

To further explore the differences between simulation and experiment in the growth 
rate of the layers, we have considered fluids made of more complicated molecules, by the 
simple device of combining Lennard-Jones atoms into chains. If the above simulations axe 
repeated with diatomic dumbbell molecules, we see no significant difference in the spreading 
behavior - at large enough solid-liquid attraction the spreading is again in the form of sharp 
layers, whose growth rate is still R 2 ~ log t. In this case however, the liquid’s volatility is 
considerably reduced so that the spreading is “dry” , without an adsorbed vapor layer on the 
substrate. For longer chains, however, the results do change in the direction of experiment. In 
work in progress [6], we have studied the spreading of chains of length 8; note that because 
of the Lennard-Jones attraction between non-adjacent atoms, the molecules are compact 
rather than elongated, and not at all entangled. Our preliminary results for this system give 
terraced spreading with a growth rate R 2 ~ t 0 8 . At the extreme, other workers [9] have 
conducted MD simulations in which the substrate is a plane whose potential varies only 
with the normal distance, which can be thought of as the limit of molecule (or monomer) 
size much larger than the solid lattice spacing. Here one finds diffusive behavior, R 2 ~ t, in 
agreement with laboratory observation. We are now studying the case of length-16 chains, 
to elucidate the crossover behavior. 

An intuitive picture of the results is as follows. The strong solid-liquid interaction draws 
fluid molecules to the solid, and the free energy is minimized by an ordered structure with 
well-defined layers, in which as many fluid atoms as possible are in positions favored by 
the periodic lattice potential. The temperature is still above melting, so the liquid atoms 
can diffuse between and within the layers, and layer growth involves vacancy creation and 
filling. If the fluid molecules are small, it is easy for most or all of their constitutent atoms 
to find the preferred sites, and the growth is slow and sub- diffusive, in analogy to diffusion 
in the presence of traps. When the molecule size increases, steric and thermal effects make 
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it difficult for a significant fraction of the atoms to simultaneously be in preferred sites, so 
that the pinning effect is less effective and the growth is faster. 

CONCLUSIONS 

In summary, we have carried out systematic molecular dynamics studies of drop spreading 
on a solid surface. We observe that fairly modest variations in the strength of the solid-liquid 
attraction potential lead to qualitatively different wetting regimes, each with a laboratory 
counterpart. We have focused on the novel terraced spreading case, where our results give 
the internal dynamical structure for this spreading regime, and where we have shown that 
the rate of spreading depends on the size and nature of the fluid molecules. 

More generally, we have indicated how molecular simulation can complement theory and 
experiment, by providing otherwise unavailable small-scale information about the internal 
dynamics of a system, by providing “clean” results free of theoretical bias of laboratory 
contamination, and by allowing one to systematically vary significant parameters in a many- 
body problem and compute the consequences. In other related work in progress, we are 
exploring the relation between the presence of structural and chemical surface heterogeneity 
and contact line hysteresis. Here again, we can control the substrate’s properties in com- 
plete detail, vary the fluid both in molecular size and in its Newtonian vs. non-Newtonian 
behavior, and observe the resulting dynamics at high resolution. We hope to help clarify 
these important issues in wetting as well. 
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1. The final states of spreading observed 
for (a) a = 1.1 at t = 635r, (b) a = 1.2 
at t = 885 r, and (c) a = 1.4 at t = 
635r. The points represent the centers 
of molecules in three-dimensions, pro- 
jected onto the front plane. 
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3. Dependance of drop radius on time 
for a = 1.4, for the first (squares) and 
the second (circles) layers. 
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2. Time evolution of complete terraced 
wetting at a = 1.4. The format is as in 
Fig. 1. 
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